home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fft / c_main.c < prev    next >
Encoding:
Text File  |  2001-08-22  |  5.7 KB  |  222 lines
  1. /* fft/c_main.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. int
  21. FUNCTION(gsl_fft_complex,forward) (TYPE(gsl_complex_packed_array) data, 
  22.                    const size_t stride, 
  23.                    const size_t n,
  24.                    const TYPE(gsl_fft_complex_wavetable) * wavetable,
  25.                                    TYPE(gsl_fft_complex_workspace) * work)
  26. {
  27.   gsl_fft_direction sign = forward;
  28.   int status = FUNCTION(gsl_fft_complex,transform) (data, stride, n, 
  29.                             wavetable, work, sign);
  30.   return status;
  31. }
  32.  
  33. int
  34. FUNCTION(gsl_fft_complex,backward) (TYPE(gsl_complex_packed_array) data,
  35.                     const size_t stride, 
  36.                     const size_t n,
  37.                     const TYPE(gsl_fft_complex_wavetable) * wavetable,
  38.                                     TYPE(gsl_fft_complex_workspace) * work)
  39. {
  40.   gsl_fft_direction sign = backward;
  41.   int status = FUNCTION(gsl_fft_complex,transform) (data, stride, n, 
  42.                             wavetable, work, sign);
  43.   return status;
  44. }
  45.  
  46. int
  47. FUNCTION(gsl_fft_complex,inverse) (TYPE(gsl_complex_packed_array) data, 
  48.                    const size_t stride, 
  49.                    const size_t n,
  50.                    const TYPE(gsl_fft_complex_wavetable) * wavetable,
  51.                                    TYPE(gsl_fft_complex_workspace) * work)
  52. {
  53.   gsl_fft_direction sign = backward;
  54.   int status = FUNCTION(gsl_fft_complex,transform) (data, stride, n, 
  55.                             wavetable, work, sign);
  56.  
  57.   if (status)
  58.     {
  59.       return status;
  60.     }
  61.  
  62.   /* normalize inverse fft with 1/n */
  63.  
  64.   {
  65.     const ATOMIC norm = ONE / (ATOMIC)n;
  66.     size_t i;
  67.     for (i = 0; i < n; i++)
  68.       {
  69.     REAL(data,stride,i) *= norm;
  70.     IMAG(data,stride,i) *= norm;
  71.       }
  72.   }
  73.   return status;
  74. }
  75.  
  76. int
  77. FUNCTION(gsl_fft_complex,transform) (TYPE(gsl_complex_packed_array) data, 
  78.                      const size_t stride, 
  79.                      const size_t n,
  80.                      const TYPE(gsl_fft_complex_wavetable) * wavetable,
  81.                      TYPE(gsl_fft_complex_workspace) * work,
  82.                      const gsl_fft_direction sign)
  83. {
  84.   const size_t nf = wavetable->nf;
  85.  
  86.   size_t i;
  87.  
  88.   size_t q, product = 1;
  89.  
  90.   TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4,
  91.     *twiddle5, *twiddle6;
  92.  
  93.   size_t state = 0;
  94.  
  95.   BASE * const scratch = work->scratch;
  96.  
  97.   BASE * in = data;
  98.   size_t istride = stride;
  99.  
  100.   BASE * out = scratch;
  101.   size_t ostride = 1;
  102.  
  103.   if (n == 0)
  104.     {
  105.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  106.     }
  107.  
  108.   if (n == 1)
  109.     {                /* FFT of 1 data point is the identity */
  110.       return 0;
  111.     }
  112.  
  113.   if (n != wavetable->n)
  114.     {
  115.       GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
  116.     }
  117.  
  118.   if (n != work->n)
  119.     {
  120.       GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
  121.     }
  122.  
  123.   for (i = 0; i < nf; i++)
  124.     {
  125.       const size_t factor = wavetable->factor[i];
  126.       product *= factor;
  127.       q = n / product;
  128.  
  129.       if (state == 0)
  130.     {
  131.       in = data;
  132.       istride = stride;
  133.       out = scratch;
  134.       ostride = 1;
  135.       state = 1;
  136.     }
  137.       else
  138.     {
  139.       in = scratch;
  140.       istride = 1;
  141.       out = data;
  142.       ostride = stride;
  143.       state = 0;
  144.     }
  145.  
  146.       if (factor == 2)
  147.     {
  148.       twiddle1 = wavetable->twiddle[i];
  149.       FUNCTION(fft_complex,pass_2) (in, istride, out, ostride, sign, 
  150.                     product, n, twiddle1);
  151.     }
  152.       else if (factor == 3)
  153.     {
  154.       twiddle1 = wavetable->twiddle[i];
  155.       twiddle2 = twiddle1 + q;
  156.       FUNCTION(fft_complex,pass_3) (in, istride, out, ostride, sign, 
  157.                     product, n, twiddle1, twiddle2);
  158.     }
  159.       else if (factor == 4)
  160.     {
  161.       twiddle1 = wavetable->twiddle[i];
  162.       twiddle2 = twiddle1 + q;
  163.       twiddle3 = twiddle2 + q;
  164.       FUNCTION(fft_complex,pass_4) (in, istride, out, ostride, sign, 
  165.                     product, n, twiddle1, twiddle2, 
  166.                     twiddle3);
  167.     }
  168.       else if (factor == 5)
  169.     {
  170.       twiddle1 = wavetable->twiddle[i];
  171.       twiddle2 = twiddle1 + q;
  172.       twiddle3 = twiddle2 + q;
  173.       twiddle4 = twiddle3 + q;
  174.       FUNCTION(fft_complex,pass_5) (in, istride, out, ostride, sign, 
  175.                     product, n, twiddle1, twiddle2, 
  176.                     twiddle3, twiddle4);
  177.     }
  178.       else if (factor == 6)
  179.     {
  180.       twiddle1 = wavetable->twiddle[i];
  181.       twiddle2 = twiddle1 + q;
  182.       twiddle3 = twiddle2 + q;
  183.       twiddle4 = twiddle3 + q;
  184.       twiddle5 = twiddle4 + q;
  185.       FUNCTION(fft_complex,pass_6) (in, istride, out, ostride, sign, 
  186.                     product, n, twiddle1, twiddle2, 
  187.                     twiddle3, twiddle4, twiddle5);
  188.     }
  189.       else if (factor == 7)
  190.     {
  191.       twiddle1 = wavetable->twiddle[i];
  192.       twiddle2 = twiddle1 + q;
  193.       twiddle3 = twiddle2 + q;
  194.       twiddle4 = twiddle3 + q;
  195.       twiddle5 = twiddle4 + q;
  196.       twiddle6 = twiddle5 + q;
  197.       FUNCTION(fft_complex,pass_7) (in, istride, out, ostride, sign, 
  198.                     product, n, twiddle1, twiddle2, 
  199.                     twiddle3, twiddle4, twiddle5, 
  200.                     twiddle6);
  201.     }
  202.       else
  203.     {
  204.       twiddle1 = wavetable->twiddle[i];
  205.       FUNCTION(fft_complex,pass_n) (in, istride, out, ostride, sign, 
  206.                     factor, product, n, twiddle1);
  207.     }
  208.     }
  209.  
  210.   if (state == 1)        /* copy results back from scratch to data */
  211.     {
  212.       for (i = 0; i < n; i++)
  213.     {
  214.       REAL(data,stride,i) = REAL(scratch,1,i) ;
  215.       IMAG(data,stride,i) = IMAG(scratch,1,i) ;
  216.     }
  217.     }
  218.  
  219.   return 0;
  220.  
  221. }
  222.